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A  computational  approach  capable  of  modeling  homogeneous  condensation  in  non-equilibrium  environments 
is  presented.  The  approach  is  based  on  the  direct  simulation  Monte  Carlo  (DSMC)  method,  extended  as 
appropriate  to  include  the  most  important  processes  of  cluster  nucleation  and  evolution  at  the  microscopic 
level.  The  approach  uses  a  recombination-reaction  energy-dependent  mechanism  of  the  DSMC  method  for  the 
characterization  of  dimer  formation,  and  the  RRK  model  for  the  cluster  evaporation.  Three-step  testing  and 
validation  of  the  model  is  conducted  by  (i)  comparison  of  clusterization  rates  in  an  equilibrium  heat  bath  with 
theoretical  predictions  for  argon  and  water  vapor  and  adjustment  of  the  model  parameters,  (ii)  comparison 
of  the  non-equilibrium  argon  cluster  size  distributions  with  experimental  data,  and  (iii)  comparison  of  the 
non-equilibrium  water  cluster  size  distributions  with  experimental  measurements.  Reasonable  agreement  was 
observed  for  all  three  parts  of  the  validation. 


I.  INTRODUCTION 

The  phenomenon  of  homogeneous  condensation  has 
regained  significant  interest  recently,  partly  because  it 
is  directly  related  to  the  nanofabrication^’^  and  partly 
because  of  the  advances  in  the  experimental®  and 
theoretical'^  methods.  Unlike  other  first-order  phase 
transitions,  the  gas-to-liquid  transition  often  occurs  in  a 
non-equilibrium  environment,  e.g.,  in  a  rapid  gas  expan¬ 
sion.  The  transient  events,  characteristic  of  the  conden¬ 
sation  in  a  non-equilibrium  environment,  are  still  difficult 
for  direct  experimental  observation,  making  validation  of 
both  empirical  correlations  and  theoretical  predictions 
challenging.  A  computational  method  capable  of  repro¬ 
ducing  diverse  experimental  condensation  scenarios,  in¬ 
cluding  both  equilibrium  and  non-equilibrium  condensa¬ 
tions,  can  provide  a  viable  alternative  to  these  methods. 

Two  different  approaches  for  modeling  the  condensa¬ 
tion  in  such  non-equilibrium  environment  as  rapidly  ex¬ 
panding  plumes  have  been  reported  in  the  literature.  The 
first  approach,  known  as  the  classical  approach,  takes  its 
starting  point  from  the  classical  nucleation  theory  (CNT) 
which  is  based  on  equilibrium  thermodynamics.®’®  The 
second  one,  known  as  the  kinetic  approach,  treats  nucle¬ 
ation  as  the  process  of  kinetic  chemical  aggregation.^ 

The  classical  approach  considers  the  energy  of  cluster 
formation  from  the  vapor  state.  Assuming  unimolecular 
reactions  of  cluster  growth  and  decay,  CNT  calculates 
the  corresponding  condensation  and  evaporation  rates 
using  the  Gibbs  distributions  and  the  principle  of  de¬ 
tailed  balance^^^®.  The  nucleation  rate  is  then  calculated 
assuming  a  steady  state  condition.^®  Although  a  rapidly 
expanding  supersonic  plume  is  quite  different  from  the 
isothermal,  ideal  gas  environment  assnmed  by  CNT,  the 
classical  predictions  were  found  to  be  qualitatively  cor¬ 


Distribution  A:  Approved  for  public 


rect  for  the  modeling  of  cluster  formation  in  supersonic 
jets.®^“^®  There  are  many  examples,  however,  when  CNT- 
based  results  cannot  be  fitted  to  experimental  data.  The 
CNT-based  prediction  of  the  cluster  size  distributions®^ 
significantly  deviated  from  experimental  data.®®  This  is 
consistent  with  the  work  of  Ref.  19,  where  it  was  also 
found  that  the  correct  prediction  of  the  cluster  size  distri¬ 
bution  along  with  the  internal  and  translational  energy 
distributions  is  beyond  the  area  of  applicability  of  the 
classical  approach. 

The  reasons  for  the  discrepancy  between  the  CNT- 
based  distributions  and  experimental  data  are  both  due 
to  problems  inherent  in  CNT  and  the  flow  conditions 
of  expanding  plumes.  The  former  include  the  ambigu¬ 
ous  definition  of  the  surface  energy  of  small  clusters,®® 
the  neglect  of  the  rotational  and  translational  degrees 
of  freedom  of  freshly  nucleated  clusters,^®  and  the  un¬ 
realistic  description  of  vapor-cluster  and  cluster-cluster 
interactions.®®  The  latter  are  related  to  the  main  as- 
snmptions  underlying  the  derivation  of  the  nucleation 
rate,  which  may  be  violated  in  rapidly  expanding  super¬ 
sonic  flows. ^®  The  transient  time  needed  for  a  system  to 
reach  steady  state  in  terms  of  the  unimolecular  cluster 
reactions  may  be  such  that  the  jet  macroparameters  will 
significantly  change  during  that  time.  Moreover,  many 
theoretical  and  experimental  results^^^^"®  suggest  that  lo¬ 
cal  thermal  equilibrium  does  not  exist  in  an  expanding 
supersonic  jet.  Thus  the  process  of  cluster  formation  is 
not  likely  to  be  isothermal.  There  have  been  recent  ad¬ 
vances  in  CNT  mainly  aimed  at  achieving  a  more  realistic 
model  for  condensation  and  evaporation  rates, 
but  other  principal  deficiencies  of  CNT  and  its  applica¬ 
tion  to  the  non-equilibrium  environment  still  have  yet  to 
be  addressed.  Note  that  despite  these  deficiencies,  the 
steady  state  CNT  nucleation  rate  with  a  correction®®  is 
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used  in  such  commercial  fluid  dynamics  codes  as  GASP 
or  CFD-particle. 

Unlike  CNT,  the  kinetic  approach  does  not  assume 
local  thermodynamic  equilibrium.  Instead,  a  micro¬ 
scopic  view  of  the  interactions  of  monomers  and  clus¬ 
ters  is  established  either  analytically  via  a  mathematical 
model,  e.g.,  by  the  Smoluchowski  equations  where  the 
interaction  between  particles  is  modeled  by  the  reaction 
rates, or  in  computer  simulations,  e.g.,  in  molecu¬ 
lar  dynamics  (MD)  calculations  where  the  interaction  is 
modeled  by  an  interaction  potential. It  can  be  shown 
that  the  application  of  the  Smoluchowski  equations  to  the 
modeling  of  clustering  in  supersonic  jets  is  computation¬ 
ally  unfeasible.  Even  though  MD  simulations  might  seem 
attractive  (since  no  information  besides  the  interaction 
potential  is  needed  to  perform  the  calculations),  they  are 
computationally  limited  to  a  system  size  of  about  one 
million  gas  particles  and  a  time  scale  of  a  few  nanosec¬ 
onds  at  most.  In  real  plumes,  such  as  thruster  nozzles  or 
ablation  jets,  the  number  of  gas  particles  and  the  expan¬ 
sion  time  is  greater  by  many  orders  of  magnitude  (see, 
for  example,  Therefore,  the  MD  technique  can¬ 

not  be  directly  applied  to  the  simulation  of  even  small 
laboratory-sized  supersonic  jets. 

A  promising  direction  in  modeling  the  coupled  con¬ 
densation  flow  is  the  use  of  a  kinetic  particle  simulation 
method,  direct  simulation  Monte  Carlo  (DSMC),^^  which 
is  applicable  in  a  wide  range  of  flow  regimes  from  free 
molecular  to  near  continuum.  It  is  a  statistical  approach 
for  solving  the  spatially  nonuniform  master  Leontovich 
equation  for  the  N-particle  distribution  function  and,  in 
the  limit  of  a  large  N,  it  represents  an  accurate  solution 
of  the  Boltzmann  equation.  The  advantage  of  DSMC  as 
compared  to  other  methods  is  that  complicated  cluster- 
cluster  and  cluster-monomer  interactions  including  the 
multi-body  reactions  of  cluster  nucleation  can  be  seam¬ 
lessly  incorporated. 

The  DSMC  method  has  been  used  to  study  the  pro¬ 
cess  of  cluster  formation  and  evolution  for  a  number  of 
years. However,  the  gas  flow  in  the  earlier  studies 
was  uniform,  the  considered  cluster  size  range  was  very 
narrow  (up  to  25  monomers  in  a  cluster)  and  the  exam¬ 
ined  reaction  types  were  unrealistically  limited  to  elastic 
collisions,  cluster  and  monomer  sticking  to  clusters,  and 
evaporation  of  monomers  from  clusters. 

More  recently,  the  DSMC  method  has  been  exten¬ 
sively  and  successfully  applied  to  modeling  the  processes 
of  cluster  formation  and  evolution  in  supersonic  jets  by 
Levin  et  al  (see,  for  example.  Refs.  16,38,39).  The  model 
initially  was  based  on  the  classical  nucleation  theory, 
with  the  new  clusters  being  formed  at  the  critical  size. 
Further  work  of  these  authors^®  extended  the  kinetic 
dimer  formation  approach  of  Ref.  41,  who  assumed  that 
a  ternary  collision  always  results  in  a  dimer  formation,  to 
include  molecular  dynamic  (MD)  simulations  for  obtain¬ 
ing  information  on  the  probability  of  dimer  formation  in 
such  ternary  collisions  The  work"^^  used  a  temperature- 
dependent  probability  of  formation  of  argon  dimers. 


Although  the  use  of  MD  has  a  number  of  advantages, 
due  to  its  inherent  limitations,  the  obtained  probability 
cannot  be  unambiguously  related  to  such  characteristics 
of  the  ternary  collision  as  the  internal  energy  of  the  col¬ 
lision  complex  and  the  kinetic  energy  of  the  impinging 
monomer.  Another  possible  limitation  of  the  above  work 
is  that  the  cell  quantities  (e.g.,  temperature)  rather  than 
individual  particle  characteristics  are  used  in  the  micro¬ 
scopic  models  of  nucleation  and  evaporation. 

In  the  present  paper,  the  DSMC  approach  for  modeling 
of  homogeneous  nucleation  in  rapidly  expanding  plumes 
is  extended  to  include  a  number  of  new  features.  Firstly, 
all  microscopic  reactions  are  modeled  using  the  charac¬ 
teristics  of  individual  particles  and  not  those  averaged 
over  DSMC  cells.  Secondly,  a  truly  kinetic  RRK  modeR^ 
is  implemented  to  characterize  the  cluster  evaporation 
rates.  Then,  an  energy  dependent  collision  procedure 
similar  to  the  recombination  reaction  model  of  Ref.  45 
is  used  for  the  collision  complex  formation.  An  empir¬ 
ical  parameter  is  used  for  the  inelastic  collision  number 
in  the  cluster-monomer  collisions.  For  dimers,  this  pa¬ 
rameter  was  calibrated  through  the  comparison  of  the 
computed  nucleation  rates  and  equilibrium  constants  in 
thermal  bath  with  available  theoretical  and  experimental 
data  for  argon  and  water.  Additional  validation  analy¬ 
sis  is  conducted  through  direct  comparison  with  reliable 
measurements  of  cluster  size  distributions  in  argon^^  and 
water^®. 


II.  COMPUTATIONAL  MODEL  OF  HOMOGENEOUS 
NUCLEATION 

In  the  present  non-equilibrium  model  of  homogeneous 
condensation  formulated  for  the  DSMC  method,  all  of 
the  most  important  processes  of  cluster  nucleation  and 
evolution  are  considered  at  the  microscopic  level.  First 
principles  of  the  kinetic  theory  are  used  to  define  the 
main  processes  of  homogeneous  condensation,  where  all 
collision,  nucleation,  and  evaporation  events  depend  on 
instantaneous  energies  of  colliding  partners,  and  not  cell 
temperature  or  other  macroscopic  quantities.  The  pro¬ 
cesses  that  are  included  in  the  model  and  described  in  de¬ 
tail  below  are  (i)  formation  of  collision  complexes  through 
the  binary  collisions  of  cluster-forming  monomer  species, 
(ii)  creation  of  dimers  through  the  collision  stabilization 
of  collision  complexes,  (iii)  elastic  monomer-cluster  colli¬ 
sions  that  change  the  translational  and  internal  energies 
of  colliding  particles,  (iv)  inelastic  monomer- cluster  colli¬ 
sions  that  result  in  monomer  sticking,  (v)  cluster-cluster 
coalescence,  (vi)  evaporation  of  monomers  from  clusters. 
The  details  on  each  of  these  processes  are  given  below. 


A.  Collision  complex  formation 

One  of  the  important  assumptions  of  the  present  model 
is  that  all  pairs  of  colliding  particles  create  collision  com- 
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plexes.  A  collision  complex  is  a  pair  of  monomers  that 
have  collided,  and  may  have  the  conditions  necessary  to 
form  a  dimer  if  struck  by  a  third  particle  during  its  life¬ 
time.  The  collision  complex  lifetime,  ti,  is  assumed  to 
be  dependent  on  the  type  of  monomers  and  their  relative 
collision  velocity,  with  the  functional  dependence  given 
by  the  well  known  Bunker’s  expression^® 

ti  =  1.5aofj,^e^E~i ,  (1) 

where  (Jq  and  Cq  are  the  potential  depth  and  separation 
distance  parameters  of  the  Lennard- Jones  (LJ)  potential, 
fjL  is  the  reduced  mass  of  the  colliding  particles,  and  E  is 
their  relative  translational  energy.  The  values  of  cro  and 
Co  used  in  this  work  are  3.2  x  10“^®  m  and  7.94  x  10“®^  J 
for  water,  and  3.405  x  10“^®  m  and  1.654  x  10“^^  J  for 
argon.  The  argon  LJ  parameters  are  well  known,  but 
H2O  values  are  roughly  estimated  by  reproducing  the 
known  viscosity  of  water  vapor  at  273  K  and  for  a  small 
range  of  temperatures  using  a  LJ  potential  model. 

The  process  of  interaction  of  collision  complexes  with 
surrounding  gas  particles  is  modeled  using  the  majorant 
frequency  scheme®®  with  the  assumption  that  the  colli¬ 
sion  complex  -  third  particle  interactions  are  governed 
by  the  Variable  Hard  Sphere  (VHS)  interaction  model.®® 
The  VHS  viscosity-temperature  exponent  w  of  a  colli¬ 
sion  complex  was  assumed  to  be  that  of  the  comprising 
monomers.  For  argon  atoms,  the  VHS  parameters  taken 
from  Ref.  35  were  assumed,  dref  =  4.17  A  and  lu  =  0.81. 
For  water  molecules,  dref  =  6.2  A  and  co  =  1  were  as¬ 
sumed,  based  on  reproducing  the  known  viscosity  of  wa¬ 
ter  vapor  at  273  K  and  for  a  small  range  of  temperatures 
using  the  VHS  potential  model.  Note  that  the  VHS  col¬ 
lision  diameters  are  used  to  compute  collision  frequency 
in  the  cell,  while  the  LJ  diameters  are  only  used  for  com¬ 
puting  the  diameter  and  lifetime  of  a  collision  complex  in 
order  to  determine  the  probability  of  a  three-body  colli¬ 
sion. 

Generally,  the  probability  P{t)  that  a  collision  com¬ 
plex  will  collide  with  a  third  particle  during  an  arbitrary 
time  r  is 


dPjr) 

dr 


=  v{1-P{t)) 


(2) 


where  v  is  the  collision  frequency  of  the  collision  complex 
with  third  particles.  For  a  mixture  of  Ng  gas  species,  v 
is  expressed  as 


Ns 

V  =  '^ni{ag),  (3) 

i=l 

where  a  is  the  corresponding  total  collision  cross-section, 
g  is  the  relative  collision  velocity  between  particles,  and 
brackets  denote  averaging  over  g.  Obviously,  the  proba¬ 
bility  that  the  collision  complex  will  move  freely  and  not 
collide  during  r  is  given  by 


This  expression  represents  the  probability  that  no  dimer 
will  be  formed  during  t.  Using  the  inverse  transform, 

Tc  =  -  h\{R)/v,  (5) 

where  R  is  a  random  number  uniformly  distributed  be¬ 
tween  0  and  1.  The  majorant  frequency  algorithm  for 
each  pair  of  colliding  monomers  is  therefore  as  follows. 

1.  Calculate  ti  and  Vmajorant  and  set  to  =  0,  z  =  1. 
Here,  Vmajorant  =  namaxgmax  is  the  majorant  collision 
frequency,  and  amax  and  gmax  are  maximum  cross  section 
and  relative  collision  velocity,  respectively. 

2.  Calculate 


,  ,  HR) 

—  ^i—1 

^majorant 


(6) 


3.  If  ti  >  ti,  go  to  the  the  next  pair  of  monomers. 

4.  If  ti  <  ti,  then  a  physical  collision  occurs  with  a 
probability  Pc  =  '^cFmaxgmax-  With  a  probability  1  - 
Pc,  i  =  i  +  \  and  the  algorithm  returns  to  step  2. 

In  this  algorithm,  for  consistency  with  the  collision 
complex  lifetime  determination,  an  expression  for  the  di¬ 
ameter  d  of  the  collision  complex  recommended  in  Ref.  49 
was  used. 


d  =  3Hoi€o/E)e,  (7) 


B.  Dimer  stabilization 


If  there  is  a  physical  collision  between  a  collision  com¬ 
plex  and  a  third  particle,  there  is  a  possibility  of  forming 
a  stable  dimer  as  a  result  of  such  a  collision.  Generally, 
the  probability  of  the  formation  of  a  stable  dimer  (dimer 
stabilization)  depends  on  the  colliding  species  and  the  en¬ 
ergies  -  both  translational  and  internal  -  of  the  colliding 
particles.  In  this  work,  constant  stabilization  probabili¬ 
ties  of  0.25  for  Ar®®  and  0.7  for  H2O  were  assumed,  which 
seem  reasonable  for  the  range  of  temperatures  under  con¬ 
sideration. 

Dimer  creation  through  the  collisional  stabilization  of 
collision  complexes  is  modeled  as  a  two-step  process,  L-l- 
M  {LM),  {LM)  +K  ^  LM  +  K.  Here,  L  and  M  are 
monomers,  (LM)  is  the  collision  complex,  and  K  is  the 
third  particle.  The  algorithm  of  this  process  is  described 
below. 

1.  Velocities  of  the  collision  complex  are  calculated  from 
the  momentum  conservation  as 


^{LM) 


niLVL  +  mMVM 
'm(LM) 


(8) 


2.  The  internal  energy  of  the  collision  complex  is  calcu¬ 
lated  from  energy  conservation, 

Elf^)=EH  +  EH  +  ^E,  (9) 


Pfreei^)  —  ^ 


(4)  where  AE  =  0.5(TO(iM)^’^LM)  “  '^-kvI  -  mLvh). 
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3.  The  total  energy  of  the  collision  complex  -  third  par¬ 
ticle  pair  is  increased  by  the  evaporation  (or  binding) 
energy  Eevap, 

Tptotal  _  rprel  ,  rpint  ,  rpint  ,  tti 

^  —  ^{LM)-K  ^(LM)  -r -^evap’ 

Here,  is  the  relative  translation  energy  of  the 

(LM)  —  K  collision.  Evaporation  energy  is  a  function 
of  cluster  size,  and  the  values  used  for  Ar  and  H2O  are 
given  in  the  following  sections.  Note  that  E"^*  may  be 
used  instead  of  since  the  vibrational  mode  of  third 
particles  will  barely  be  excited  at  the  low  temperatures 
at  which  homogeneous  condensation  usually  occurs. 

4.  New  energies  sampled  us¬ 

ing  the  Larsen-Borgnakke  scheme^^  extended  to  multiple 
energy  modes. 


C.  Elastic  and  inelastic  reflective  collisions  of  monomers 
and  clusters 

The  collisions  between  monomers  and  clusters  are  one 
of  the  key  processes  that  determine  the  nucleation  rate. 
The  reason  for  this  is  the  strong  dependence  of  the  evap¬ 
oration  rate  on  the  cluster  internal  energy.  Since  the 
monomers  are  dominant  in  the  flows  considered  here,  the 
cluster  internal  energy  is  mostly  governed  by  its  relax¬ 
ation  through  cluster-monomer  collisions.  In  this  work, 
a  hard  sphere  model  is  assumed  for  cluster-monomer  col¬ 
lisions,  with  the  cluster  diameter  determined  from  Eqn.  7 
for  dimers,  and  for  larger  clusters  from  an  empirical  cor¬ 
relation  used  extensively  in  the  past  (see,  for  example, 
Ref.  16), 

d  =  2-{A-i^  +B),  (11) 

where  A  and  B  are  species-dependent  constants,  and  i  is 
the  number  of  monomers  in  the  cluster.  In  this  work,  the 
values  of  A  and  B  were  2.3  x  10“^°  m  and  3.4  x  10“^°  m 
for  argon, and  1.9  x  10“^°  m  and  2.4  x  10“^°  m  for 
water. 

For  the  energy  transfer  between  the  relative  transla¬ 
tional  and  internal  modes  of  the  cluster  and  monomer, 
the  Larsen-Borgnakke  model®^  is  used,  and  a  parameter 
Z  is  introduced  which  has  a  meaning  of  the  internal  en¬ 
ergy  relaxation  number.  The  energy  transfer  between  all 
energy  modes  of  the  cluster-monomer  pair  occurs  with  a 
probability  Z~^,  and  an  elastic  collision  with  no  inter¬ 
nal  energy  exchange  occurs  with  the  additional  proba¬ 
bility  1  —  Z~^.  For  argon,  temperature  dependent  values 
of  Z(T)  were  used,  obtained  through  the  linear  interpo¬ 
lation  between  the  values  given  in  Table  I.  Similar  to 
temperature  dependent  collision  numbers  for  rotational 
and  vibrational  relaxation  of  molecules  widely  used  in  the 
DSMC  method,  T  was  the  cell-based  translational  tem¬ 
perature.  Such  a  temperature  dependence  allows  good 
agreement  of  the  DSMC  rates  for  dimer  nucleation  and 
dissociation  with  rates  available  in  the  literature.  For  wa¬ 
ter  condensation,  a  constant  value  of  Z  =  10  was  used. 


TABLE  I.  Inelastic  collision  number  as  function  of  tempera¬ 
ture  for  argon. 


T,  K 

0.0 

100.0 

200.0 

300.0 

400.0 

500.0 

0.25 

0.13 

0.08 

0.06 

0.046 

0.04 

TABLE  II.  Sticking  probability  as  function  of  cluster  size  for 
argon. 


Cluster  size 

Probability 

2 

0.06 

3 

0.075 

4 

0.1 

5 

0.16 

6 

0.3 

7 

0.5 

8 

0.67 

9 

0.75 

10 

0.8 

11 

0.833 

12 

0.857 

13 

0.876 

14 

0.891 

15 

0.9 

which  corresponds  to  that  for  rotational  and  vibrational 
relaxation  of  water  molecules  at  room  temperature  ob¬ 
tained  in  molecular  dynamics  studies. For  argon,  the 
values  of  Z  given  in  Table  I  are  used,  based  on  the  fitting 
described  in  Section  III. 


D.  Sticking  collisions  of  monomers  and  clusters 


When  a  monomer  collides  with  a  cluster,  sticking  of 
the  monomer  to  the  cluster  surface  is  possible,  in  addi¬ 
tion  to  a  reflective  collision  described  in  the  previous  sec¬ 
tion.  For  small  clusters,  monomer  sticking  is  the  main 
process  that  governs  the  evolution  of  the  droplet  size 
distribution."^®  For  water  molecules,  an  empirical  depen¬ 
dence  of  the  sticking  probability  on  the  species  radius  and 
mass,  given  in  Ref.  56,  is  used.  This  dependence  reduces 
to 


e 


dl  f  rn„  \  " 
{dn  +  di)^  \mn  +  mi) 


(12) 


where  indices  n  and  I  refer  to  the  cluster  and  monomer, 
respectively.  Note  that  for  low  n  this  expression  agrees 
well  with  the  molecular  dynamics  results  of  Ref.  57.  For 
argon,  the  sticking  probability  of  monomers  on  clusters 
given  in  Table  II  is  used,  based  on  the  results  of  molecular 
dynamics  simulations  of  Ref.  38. 

The  following  algorithm  is  used  to  model  sticking  of  a 
molecule  to  a  cluster  L  to  form  a  larger  cluster  K. 

I.  Calculate  velocities  of  K  as 

rriLVL  +  mMVM 

Vk  =  -  (13) 

nik 
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2.  After  velocities  vk  are  assigned,  the  new  internal  en¬ 
ergy  is  calculated 

+  E^^*  +  Q  +  l\E  (14) 

Here, 

l\E  =  -0.5{mKVK  +  rriLvl  +  itimVm) 

and  Q  =  E^^ap  is  the  evaporation  energy  of  one  monomer 
off  cluster  K. 

Similar  to  monomer-cluster  collisions,  the  outcome  for 
cluster-cluster  collisions  is  assumed  to  be  either  coales¬ 
cence  or  elastic  interaction.  The  probability  of  stick¬ 
ing  was  assumed  to  be  unity  both  for  argon  and  water 
in  cluster-cluster  collisions.  The  algorithm  for  cluster- 
cluster  sticking  collisions  is  similar  to  the  monomer- 
cluster  collision,  with  the  exception  of  Q  =  —Qk  +  Ql  + 
Qmi  where  Qi  is  the  energy  of  vaporization  of  cluster  i. 


E.  Evaporation  rate  and  algorithm 


F.  Energy  redistribution  in  evaporation 

The  following  energy  redistribution  scheme  is  used  to 
model  the  evaporation  of  a  monomer  M  off  a  cluster  K, 
with  a  daughter  cluster  L  formed. 

1.  Decrease  the  internal  energy  of  K  by  the  evaporation 
energy  Eevap 

T^int  _  -ppint  /'17^ 

~  ~  -^evap 

2.  Split  cluster  internal  energy  between  the  relative 
translational  energy  Erei  and  the  internal  energy  of  the 
M  -  L  pair  (the  sum  of  the  daughter  cluster  internal 
energy  and  the  monomer  internal  energy) ,  using  the 
Larsen-Borgnakke  procedure. 

3.  Split  energy  Ei„t  between  the  daughter  cluster  inter¬ 
nal  energy  and  the  monomer  internal  energy  using  the 
Larsen-Borgnakke  procedure. 

4.  Calculate  new  daughter  cluster  and  monomer  veloci¬ 
ties  using  Erei,  keeping  in  mind  that  the  center  of  mass 
of  the  new  pair  has  the  same  velocity  as  the  velocity  of 
the  original  cluster. 


Following  Ref.  34,  RRK  theory  is  used  to  model  the 
evaporation  process.  The  evaporation  rate  kg  is  calcu¬ 
lated  as 


G.  Evaporation  energy  and  the  number  of  internal 
degrees  of  freedom 


ke 


vNs 


-  E 
Eint 


evap 


3n-7 


(15) 


Here,  n  is  the  number  of  monomers  in  the  cluster,  v  is  the 
vibration  frequency,  Ng  is  the  number  of  surface  atoms, 
Eevap  is  the  evaporation  energy,  and  Eint  is  the  cluster  in¬ 
ternal  energy.  For  dimers,  the  exponent  3n  — 7  is  replaced 
with  1.  The  number  of  surface  atoms  Ng  is  n  for  n  <  5, 
n  —  1  for  4  <  n  <  7,  and  (367r)^/^(n^/^  —  1)^  for  n  >  6. 
The  vibration  frequency  was  taken  to  be  2.68  x  10^^  s“^ 
for  water  clusters,^®  and  10^^  s“^  for  argon  clusters.®® 
Since 


^  =  -keNg,  (16) 

at 

the  time  r  to  the  next  evaporation  event  for  a  given  clus¬ 
ter  may  be  sampled  from  r  =  —  In  ^ .  The  algorithm 
used  to  model  the  evaporation  process  for  a  given  cluster 
over  a  simulation  timestep  At  is  as  follows. 

1.  Set  tiocai  —  d 

2.  Calculate  kg 

3.  Change  tiogai  —  1/oca/— i  Iti(^R^ / kg 

4.  If  tiocai  >  At,  exit. 

5.  Evaporate  one  monomer  (see  below). 

6.  If  the  remaining  cluster  is  a  monomer,  exit.  Else,  go 
to  Step  3. 

Note  that  in  a  limiting  case  of  large  clusters  (cluster 
size  tends  to  infinity),  Eqn.  (15)  reduces  to  an  exponen¬ 
tial  expression  similar  to  that  recently  obtained  using  the 
unimolecular  evaporation  theory.®® 


It  is  important  to  use  reasonable  values  for  the  evapo¬ 
ration  energy  of  a  monomer  off  a  cluster,  as  well  as  the 
number  of  cluster  internal  degrees  of  freedom,  as  a  func¬ 
tion  of  cluster  size.  The  reason  for  this  is  that  all  energy 
redistribution  processes  that  involve  clusters  depend  di¬ 
rectly  on  the  number  of  internal  degrees  of  freedom.  The 
evaporation  energy  of  a  monomer  from  an  n-cluster  is 
given  by  E},{n)  —  Eb{n  —  1),  where  Ej,  is  the  binding  (en¬ 
ergy  to  separate  all  n  monomers  in  an  n-cluster)  energy  of 
the  n-cluster.  The  Et  values  for  water  clusters  (Ref.  60) 
were  taken  from  Ref.  61  for  n  =  2  —  8  (and  subtracting 
the  zero-point  energies,  taken  from  Ref.  62),  and  using 
smoothed  values  from  Refs.  63,64  for  n  =  9  —  13.  Evap¬ 
oration  energies  were  taken  from  Ref.  65  for  argon  clus¬ 
ters.  The  number  of  cluster  internal  degrees  of  freedom 
is  calculated  from  the  expression  for  the  average  internal 
energy  (E)  of  a  cluster  of  a  size  n 


as 


{E) 


eint  o 

^kT  =  nCvT  -  -kT 
2  2 


int 


-  3 


where  Cy  is  the  cluster  heat  capacity.  For  water,  the 
values  of  Cy  where  adapted  from  Refs.  66,67.  For  water 
dimers,  the  value  of  Cy  was  taken  from  Ref.  67  at  200  K, 
and  the  number  of  degrees  of  freedom  was  assumed  to  be 
a  function  of  internal  energy  (or  effective  temperature), 
decreasing  linearly  from  its  listed  value  at  200  K  to  3  at 
0  K.  The  heat  capacity  for  n  =  8  and  n  =  10  was  taken 
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TABLE  III.  The  water  and  argon  cluster  heat  capacities  and 
evaporation  energies  per  monomer. 


Size 

H2O  Eevapt  J 

H2O  Cy,  J/K 

Ar  Eevapi  J 

Ar  Cy,  J/K 

2 

2.455E-20 

3.726E-23 

1.98E-21 

2.41E-23 

3 

5.352E-20 

4.830E-23 

3.96E-21 

2.76E-23 

4 

6.325E-20 

5.748E-23 

5.94E-21 

3.10E-23 

5 

4.726E-20 

6.555E-23 

6.08E-21 

3.31E-23 

6 

4.587E-20 

7.245E-23 

6.89E-21 

3.44E-23 

7 

5.560E-20 

7.817E-23 

7.49E-21 

3.54E-23 

8 

6.950E-20 

8.277E-23 

6.36E-21 

3.62E-23 

9 

5.560E-20 

8.648E-23 

8.29E-21 

3.68E-23 

10 

5.699E-20 

8.970E-23 

8.25E-21 

3.72E-23 

11 

5.838E-20 

9.227E-23 

8.27E-21 

3.76E-23 

12 

5.977E-20 

9.467E-23 

9.86E-21 

3.79E-23 

13 

6.116E-20 

9.654E-23 

12.2E-21 

3.82E-23 

from  Ref.  66  at  200  K.  The  values  of  Cy  were  interpolated 
between  n  =  2,  8,  10,  and  the  bulk  liquid  value.  For 
argon,  the  values  of  Cy  were  assumed  to  approximate  an 
expression  =  2(3n  —  rj)  +  e,  where  rj  =  5  and  e  =  2 
for  n  =  2  and  rj  =  6  and  e  =  3  otherwise.®®  Both  the  heat 
capacity  and  evaporation  energies  are  listed  in  Table  III. 
For  the  cluster  sizes  larger  than  given  below,  the  values 
for  the  maximum  listed  sizes  are  used. 

The  non-equilibrium  condensation  model  described 
here  was  implemented  in  the  DSMC  code  SMILE®®.  The 
validation  of  the  code  through  the  comparison  with  the¬ 
oretical  and  experimental  results  is  presented  below. 

III.  THERMAL  BATH  RELAXATION 

Inelastic  cross  sections  for  monomer-monomer  and 
monomer- cluster  collision  processes  are  necessary  for 
detailed  validation  of  a  kinetic  condensation  model. 
These  cross  sections,  generally  functions  of  the  energy 
states,  both  translational  and  internal,  of  pre-  and  post- 
collisional  particles,  are  not  available  for  most  gas  and 
temperature  conditions  of  interest.  Contrary  to  the  en¬ 
ergy  dependent  cross  section,  the  integral  temperature 
dependent  rates  for  such  collisions  at  conditions  close 
to  equilibrium  are  available  in  the  literature  for  argon 
dimers.  Therefore,  one  of  the  key  indicators  of  the  accu¬ 
racy  and  reliability  of  a  condensation  model  is  its  ability 
to  produce  realistic  rates  of  evaporation  and  nucleation 
at  eqnilibrium.  Although  matching  the  rates  generally 
does  not  guarantee  correct  behavior  in  nonequilibrium, 
it  still  is  a  necessary  condition  for  a  model  to  satisfy. 

In  this  work,  thermal  bath  relaxation  of  pure  argon  and 
pure  water  vapor  are  examined  at  different  temperature 
conditions  with  the  DSMC  method,  and  the  dimer  forma¬ 
tion  rates  for  argon  and  equilibrium  constants  for  the  for¬ 
mation  of  dimers  in  argon  and  water  are  calculated  and 
compared  to  the  published  results. In  all  thermal 
bath  results,  one  million  simulated  particles  were  used, 
and  the  run  proceeded  until  the  steady  state  is  reached, 
after  which  the  results  were  sampled  for  20  thousand 


FIG.  1.  Argon  dimer  formation  rate  as  a  function  of  gas 
temperature. 

timesteps.  The  number  density  was  5  x  10®®  molec/m® 
for  argon  and  2  x  10®®  molec/m®  for  water.  The  timestep 
of  2.5  X  10“®®  s  was  selected  so  that  the  number  of  col¬ 
lisions  per  molecule  is  much  smaller  than  unity,  and  the 
results  are  independent  on  the  timestep. 

The  computed  dimer  formation  rates  kyec  for  argon  are 
presented  in  Fig.  1  and  compared  with  the  stable  dimer 
formation  rate  of  Ref.  70,  where  they  were  calculated 
using  classical  trajectories,  and  the  following  expression 
was  proposed, 

krec  =  10.15T-O  ®®'®  exp  {-O.OOdlT}  .  (18) 

Generally,  the  present  dimer  formation  rates  are  in  good 
agreement  with  the  classical  trajectory  calculations,  with 
the  maximum  difference  approaching  20%  for  higher  tem¬ 
peratures.  Note  that  the  computed  rate  has  a  some¬ 
what  different  slope  than  that  of  Eqn.  18.  A  number 
of  factors  could  be  affecting  the  slope,  among  whieh  are 
the  energy  dependence  of  the  stabilization  probability 
and  dimer  heat  capacity,  which  were  not  included  in  the 
present  model.  The  MD  calculations  in  Ref.  40  for  ar¬ 
gon  dimer  formation  predicted  a  small  decrease  in  sta¬ 
bilization  probability  with  increasing  temperature  over  a 
limited  temperature  range. 

The  computed  equilibrium  constant,  K^q,  which  is  the 
ratio  of  the  dimer  dissociation  to  the  dimer  formation 
rate,  is  given  in  Fig.  2.  It  is  compared  to  the  theoret¬ 
ical  results  of  Ref.  71,  where  a  number  of  approximate 
classical  and  quantum  methods  are  compared  with  exact 
numerical  calculations,  and  also  experimental  results  of 
Ref.  72  and  theoretical  results  of  Ref.  73.  Note  that  while 
the  results  for  different  models  and  interaction  potentials 
were  found  to  be  widely  different  in  Ref.  71,  there  was  a 
good  agreement  between  analogous  quantum  and  classi- 
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FIG.  2.  Argon  dimer  formation  equilibrium  constant  as  a 
fnnction  of  gas  temperature. 


cal  calculations.  Figure  2  shows  that  for  the  entire  tem¬ 
perature  range  there  is  an  excellent  agreement  between 
the  present  model  and  the  theoretical  and  experimental 
values.  The  reason  for  such  a  good  agreement  is  the  ap¬ 
propriate  selection  of  the  temperature  dependence  of  the 
inelastic  collision  number  Z . 

This  parameter,  which  is  in  effect  the  inverse  probabil¬ 
ity  of  the  energy  transfer  between  the  internal  modes  of 
a  dimer  and  the  translational  modes  in  dimer-monomer 
collisions,  was  found  to  be  an  important  factor  that  in¬ 
fluences  the  magnitude  of  the  equilibrium  constant  K^^q. 
This  may  be  explained  as  follows.  The  dimers  are  formed 
after  three-body  collisions,  and  typically  have  internal 
energies  smaller  than  the  evaporation  energy  after  those 
collisions.  In  argon,  the  evaporation  energy  for  a  dimer 
is  relatively  small  compared  to  the  typical  total  col¬ 
lision  energy  for  all  temperatures  under  consideration 
{Egyap/k  ~  140  K).  That  means  that  most  of  the  dimers 
will  have  their  internal  energy  in  excess  of  the  evapora¬ 
tion  energy  just  after  one  or  two  inelastic  collisions  with 
monomers  for  thermal  bath  temperatures  greater  than 
140  K.  The  lifetime  of  the  dimers  whose  internal  energy 
is  larger  than  the  evaporation  energy  is  very  short,  on  the 
order  of  a  picosecond.  This  results  in  the  dimer-monomer 
energy  transfer  being  the  main  process  that  leads  to  quick 
dimer  dissociation.  Note  that  the  value  of  Z  has  negli¬ 
gible  impact  on  the  dimer  formation  rates,  and  only  the 
evaporation  rates  are  affected.  As  a  result,  in  the  range 
of  temperatures  considered  in  this  work,  the  equilibrium 
constant  for  argon  was  found  to  be  nearly  proportional 
to  Z~^. 

The  Z  dependence  of  the  equilibrium  constant  is 
weaker  for  water  condensation.  In  this  case,  the  evap¬ 
oration  energy  of  a  dimer  is  much  larger  than  the  trans¬ 


lational  energy  of  colliding  molecules  and  dimers  (the 
reduced  evaporation  energy  Egyap/k  «  1,800  K,  com¬ 
pared  to  gas  temperatures  on  the  order  of  300  K).  The 
high  value  of  the  evaporation  energy  results  in  longer  life¬ 
times  of  dimers,  since  many  more  collisions  are  necessary 
to  transfer  enough  energy  from  the  translational  modes 
to  the  internal  modes  of  a  dimer.  The  dependence  of  Kgq 
on  Z  is  therefore  much  weaker  for  water  than  for  argon. 
In  a  250  K  thermal  bath,  K^g  was  found  to  decrease  by 
only  about  a  factor  of  two  when  Z  decreases  from  10  to 
1. 

Comparison  of  the  equilibrium  constant  obtained  with 
the  present  model,  with  the  theoretical  results  of  Refs.  67, 
74,  where  a  flexible  potential  energy  surface  fitted  to 
spectroscopical  data  was  used,  is  shown  in  Fig.  3.  Note 
that  there  are  a  number  of  theoretical  predictions  of  the 
equilibrium  rate  of  water  dimerization,  and  they  differ 
by  at  least  a  factor  of  three  in  the  range  of  temperatures 
considered  in  this  work.  Reference®^  was  chosen  for  com¬ 
parison  as  the  most  sophisticated  and  one  of  the  most 
recent  ones.  The  results  show  that  the  calculated  equilib¬ 
rium  constant  agrees  well  with  the  theoretical  prediction 
for  temperatures  between  160  K  and  250  K,  the  range 
that  is  expected  to  be  very  important  in  terms  of  dimer 
formation  and  nucleation  in  plume  expansions.  The  cal¬ 
culated  values  for  higher  temperatures  significantly  over¬ 
predict  the  theoretical  curve.  It  is  important  to  note  that 
these  results  are  sensitive  to  the  values  of  Cy  and  Egyap- 
For  example,  a  20%  change  in  Cy  results  in  a  factor  of 
two  difference  in  Kgq.  However,  the  computations  have 
shown  that  the  change  in  the  binding  energy  and  heat  ca¬ 
pacity,  that  are  assumed  to  be  temperature  independent, 
do  not  significantly  change  the  slope  of  the  equilibrium 
constant  decreasing  with  increasing  temperature.  This 
general  trend  of  weaker  slope  may  be  related  to  a  num¬ 
ber  of  factors,  such  as  temperature  dependence  of  ac¬ 
tual  heat  capacity  and  evaporation  energy,  as  well  as  the 
approximate  after-reaction  energy  redistribution  used  in 
this  work  (recall  that  the  Larsen-Borgnakke  model  was 
used  for  the  energy  redistribution) .  Another  possible  and 
likely  reason  is  the  temperature  dependence  of  the  dimer 
stabilization  probability  for  three-body  collisions,  that 
was  not  included  in  the  present  model. 


IV.  ARGON  CLUSTER  SIZE  DISTRIBUTION  IN 
ORIFICE  EXPANSION 

The  ability  of  a  condensation  model  to  reasonably  cap¬ 
ture  theoretical  dimer  equilibrium  constant  and  dimer 
formation  rate  at  equilibrium  is  crucial  for  accurate  pre¬ 
diction  of  the  amount  of  dimers  in  a  condensing  gas 
flow.  However,  it  does  not  guarantee  correct  rates  for  nu¬ 
cleation  and  evaporation  of  clusters  larger  than  dimers. 
Therefore,  in  order  to  assess  the  model  performance,  it 
is  necessary  to  compare  cluster  size  distributions,  ob¬ 
tained  in  numerical  simulations,  with  available  experi¬ 
mental  data.  Measurements  of  the  pressure  dependence 


FIG.  3.  Water  equilibrium  constant  as  a  function  of  gas  tem¬ 
perature. 


of  size  selected  argon  clusters  present  a  good  validation 
tool  for  DSMC  based  condensation  modeling.  The  rea¬ 
sons  for  selecting  these  measurements  is  that  their  oper¬ 
ating  Knudsen  numbers,  while  rather  low,  are  still  high 
enough  for  the  DSMC  method  to  be  used,  and  their  accu¬ 
racy  is  quite  high.  The  results  are  based  on  the  method  of 
size  selection  of  the  clusters  in  a  scattering  process  with 
He/sde  Originally  designed  to  measure  reliable  fragmen¬ 
tation  probabilities  in  the  detection  process,  the  data  can 
also  be  used  to  measure  the  fractions  of  different  cluster 
sizes  as  function  of  the  pressure  of  the  expansion. 

The  original  data  were  complemented  by  new  results  up 
to  n  =  10  and  newly  evaluated  using  the  recently  ob¬ 
tained  fragmentation  probabilities  in  this  size  range. 
The  results  for  the  averaged  cluster  sizes  appeared  al¬ 
ready  in  Ref.  47. 

In  order  to  reproduce  the  experimental  setup,  an  ar¬ 
gon  expansion  out  of  a  thin  40  fxm  diameter  orifice  was 
modeled,  with  a  stagnation  and  orifice  wall  tempera¬ 
tures  of  300  K.  Two  stagnation  pressures  were  consid¬ 
ered,  1.5  atm  and  2  atm.  The  simulated  part  of  the 
plenum  was  70  /rm  in  axial  and  radial  directions,  with  the 
equilibrium  stagnation  conditions  imposed  at  the  inflow 
boundaries.  The  size  of  the  subsonic  part  was  chosen  suf¬ 
ficiently  large  to  minimize  the  impact  of  the  orifice  on  the 
inflow  boundaries.  The  plume  part  of  the  computational 
domain  was  stretched  to  700  fxm  in  the  axial  direction  to 
capture  the  terminal  argon  cluster  mole  fractions,  similar 
to  those  recorded  in  the  experiments.  Vacuum  boundary 
conditions  were  imposed  on  the  outflow  boundaries.  The 
the  axis  of  symmetry  boundary  condition  was  imposed  at 
the  lower  boundary.  Fully  diffuse  accommodation  on  the 
wall  was  assumed.  To  avoid  slow  convergence  associated 
with  the  subsonic  part  of  the  domain,  the  calculations 


included  two  steps.  First,  the  flow  was  calculated  in  the 
free  molecular  regime.  Then,  the  simulated  molecules 
from  the  first  step  were  utilized  as  the  initial  condition 
for  the  second,  high  pressure,  step.  The  numbers  of  sim¬ 
ulated  molecules  and  collision  cells  were  8  million  and  1.5 
million,  respectively. 


The  gas  temperature  field  for  the  1.5  atm  stagnation 
pressure  case  is  shown  in  Fig.  4  (left).  Only  part  of  the 
computational  domain  is  shown  to  provide  more  detail 
to  the  plume  near  field.  The  gas  density  in  the  plenum  is 
fairly  high,  which  causes  relatively  large  number  of  col¬ 
lisions  in  the  plume.  The  latter  results  in  low  freezing 
temperatures,  about  4  K  along  the  nozzle  axis.  The  flow 
starts  to  freeze  after  about  500  /xm  downstream  from  the 
orifice  exit  plane.  The  density  of  the  expanding  plume 
keeps  its  sharp  decrease,  as  shown  in  Fig.  4  (right).  The 
decreasing  number  of  collisions  in  the  flow,  and  especially 
three-body  collisions  that  govern  the  dimer  formation 
rate,  results  in  freezing  the  mole  fraction  of  argon  clusters 
at  500  fim,  which  is  also  shown  in  Fig.  4  (right).  Note 
that  a  smoothed  mole  fraction  of  all  clusters  is  given  in 
this  hgure.  In  the  plenum,  it  is  mostly  dimers;  their  num¬ 
ber  in  that  region  corresponds  to  the  equilibrium  mole 
fraction  at  a  given  temperature  and  density. 


Comparison  of  computed  and  measured  cluster  size 
distributions  for  the  two  pressure  under  consideration  is 
presented  in  Fig.  5.  The  relative  mole  fractions  are  given, 
with  the  combined  total  adding  up  to  the  unity.  For  the 
1.5  atm  case,  there  is  a  very  good  agreement  between  nu¬ 
merical  modeling  and  experiment  for  dimers  and  trimer. 
For  larger  clusters,  the  computed  values  are  lower  than 
measured.  Note  that  the  statistical  nature  of  the  DSMC 
method  makes  it  very  difficult  to  calculate  the  concen¬ 
tration  of  cluster  sizes  which  relative  mole  fraction  is 
less  than  0.01,  with  acceptable  statistical  accuracy.  Even 
though  each  computation  took  days  on  a  multi-processor 
computer,  the  values  of  the  error  bars  for  such  cluster 
sizes  are  still  on  the  order  of  the  values  of  the  correspond¬ 
ing  mole  fractions.  The  same  is  true  for  the  experimental 
errors  for  the  small  fractions.  Accounting  for  these  nu¬ 
merical  and  experimental  error  bars,  it  may  be  concluded 
that  there  is  a  satisfactory  agreement  between  the  model¬ 
ing  and  the  measurement  for  the  1.5  atm  case.  For  2  atm, 
the  mole  fractions  of  the  clusters  larger  than  trimers  are 
noticeably  higher  than  for  1.5  atm,  and  the  agreement 
between  the  DSMC  and  the  experiment  is  better.  Again, 
the  computational  and  experimental  points  agree  within 
the  limits  of  uncertainties  of  the  used  techniques  and 
models.  Note  also  that  the  computed  dimer  mole  fraction 
is  in  reasonable  agreement  with  the  Knuth’s  semiempiric 
expressions^®.  Knuth’s  non-dimensional  scaling  parame¬ 
ter,  that  defines  the  dimer  mole  fraction,  is  0.0037  and 
0.006  as  compared  to  the  DSMC  values  of  0.0053  and 
0.009  for  pressures  of  1.5  atm  and  2  atm,  respectively. 


9 


FIG.  4.  Gas  temperature  field  (left)  and  number  density  and  cluster  mole  fraction  profiles  along  the  orifice  axis  (right)  in  argon 
orifice  expansion.  1.5  atm  case. 


FIG.  5.  Terminal  argon  cluster  relative  mole  fraction  for  a  stagnation  pressure  of  1.5  atm  (left)  and  2.0  atm  (right). 


V.  WATER  CLUSTER  SIZE  DISTRIBUTION  IN 
NOZZLE  FLOW 

The  last  part  of  the  validation  and  numerical  analy¬ 
sis  of  the  presented  condensation  model  is  focused  on 
the  nucleation  and  evolution  of  small  water  clusters  in 
a  conical  nozzle.  The  study  was  prompted  by  the  avail¬ 
ability  of  high  quality  experimental  data"^®  on  terminal 
size  distribution  of  water  clusters  in  the  range  of  flow 
conditions  where  the  pressure  was  relatively  low  so  that 
the  computational  cost  of  using  the  DSMC  method  is  not 
prohibitive  (although  still  rather  high).  The  results  are 
obtained  by  doping  the  water  clusters  by  one  Na  atom, 
which  is  photoionized  close  to  the  threshold  without  frag¬ 
mentation.  The  results  are  presented  for  the  nozzle  ge¬ 
ometry  of  Ref.  48.  The  nozzle  is  a  conical  nozzle  with  a 


41°  opening  angle,  a  total  length  of  2  mm,  and  a  throat 
diameter  of  50  /xm.  The  two  smallest  stagnation  pres¬ 
sures  considered  in  Ref.  48,  1.577  bar  and  2.173  bar,  are 
used  in  this  work,  with  the  corresponding  stagnation  tem¬ 
perature  of  495  K.  A  constant  nozzle  wall  temperature 
of  495  K  was  used  to  reproduce  the  experimental  setup. 
Since  the  background  pressure  effect  in  the  experiment  is 
believed  to  be  small, the  flow  expansion  into  the  vac¬ 
uum  is  modeled. 

The  axisymmetric  capability  of  SMILE  was  used,  with 
the  total  number  of  simulated  molecules  and  collision 
cells  about  80  million  and  10  million,  respectively,  for 
a  pressure  of  1.577  bar,  and  160  million  and  20  million 
for  a  pressure  of  2.173  bar.  A  uniform  400x100  grid  was 
used  for  sampling  of  macroparameters  and  distribution 
functions.  The  Larsen-Borgnakke  model  with  constant 
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relaxation  numbers  of  10  for  both  rotational  and  vibra¬ 
tional  energies  was  used  for  energy  transfer  in  monomer- 
monomer  collisions,  and  the  reflection  of  particles  on  the 
nozzle  surface  was  assumed  to  be  fully  diffuse.  Uniform 
inflow  conditions  were  imposed  at  the  nozzle  throat,  cal¬ 
culated  from  the  isentropic  flow  relations. 

The  first  set  of  results  presented  here  shows  the  effect 
of  the  condensation  on  the  gas  flow  inside  the  nozzle  and 
in  the  plume  near  field.  The  gas  translational  tempera¬ 
ture  and  axial  velocity  is  shown  in  Fig.  6  for  two  cases, 
the  baseline  condensation  model  and  the  condensation 
turned  off.  The  results  show  that  the  condensation  prac¬ 
tically  does  not  change  the  flow  parameters  inside  the 
boundary  layer.  This  is  expected,  since  the  temperature 
in  the  boundary  layer  is  higher  than  in  the  coreflow,  and 
nucleation  near  the  surface  is  not  likely.  Near  the  cen¬ 
terline,  though,  the  condensation,  being  an  exothermic 
process,  results  in  a  significant  heat  release.  The  tem¬ 
perature  in  that  region  visibly  increases,  with  the  max¬ 
imum  change  of  over  20  K.  The  higher  temperature  for 
the  condensing  flow,  accompanied  by  fast  translational 
relaxation,  causes  an  increase  in  the  axial  velocity  in  the 
coreflow  of  about  30  m/s.  The  velocity  change  inside  the 
boundary  layer  is  negligibly  small. 

The  computations  conducted  for  a  higher  stagnation 
pressure  of  2.173  bar  show  visible  effects  of  pressure  on 
both  gas  temperature  and  number  of  clusters  in  the  flow 
(see  Fig.  7).  As  expected,  the  boundary  layer  thickness 
decreases  with  increasing  pressure.  The  temperatures  in 
the  coreflow  are  similar  only  during  the  first  100  /xm  from 
the  nozzle  throat.  Further  downstream,  mostly  due  to 
the  condensation  heat  release,  the  temperature  in  the 
coreflow  is  up  to  15  K  higher  for  the  larger  pressure  case. 
The  cluster  mole  fraction  is  also  noticeably  higher.  Gen¬ 
erally,  the  cluster  concentration  increases  rapidly  in  the 
first  100/xm  from  the  nozzle  throat,  and  reaches  its  max¬ 
imum  at  about  300/im  for  the  1.577  bar  case  and  200/im 
for  2.173  bar.  After  that,  the  evaporation  and  cluster 
coalescence  result  in  some  decrease  of  the  cluster  concen¬ 
tration.  The  maximum  mole  fraction  is  1.5%  for  the  lower 
pressure  case,  and  2.5%  for  the  higher  pressure  case.  In 
addition  to  an  increase  in  the  number  of  clusters  with 
stagnation  pressure,  the  average  cluster  size  also  slightly 
increases  with  pressure,  although  this  increase  is  much 
less  pronounced  that  that  of  the  total  number  of  clus¬ 
ters.  The  average  terminal  cluster  size  is  about  4.1  and 
4.2  for  1.577  bar  and  2.173  bar,  respectively. 

Comparison  of  the  terminal  size  distribution  obtained 
using  the  presented  model,  with  experimental  results'^®  is 
shown  in  Fig.  8.  In  the  experiments,  the  dimers  were  be¬ 
low  the  detection  threshold,  the  trimer  population  may 
be  somewhat  affected  by  that  threshold,  and  all  larger 
clusters  are  believed  to  be  recorded  without  significant 
distortions.  Exponential  prohles  of  the  size  distribution 

function,  /  oc  exp  >  that  have  the  average  cluster 

size  (n)  that  corresponds  to  the  experimental  values  of 
7  and  20,  are  also  shown  in  this  figure.  For  side-by-side 
comparison,  the  experimental  size  distribution  was  nor¬ 


malized  so  that  it  has  the  same  fraction  of  clusters  larger 
than  dimers  as  in  the  exponential  distribution  approxi¬ 
mating  it. 

The  results  show  that  there  is  a  reasonable  agreement 
between  the  numerical  modeling  and  measurements  for 
the  lower  pressure.  The  computed  average  cluster  size  for 
1.577  bar  was  5.5,  as  compared  to  the  measured  value 
of  7±2.  This  value  is  based  on  a  reevaluation  of  the 
data  by  correctly  subtracting  the  background.  There  is 
a  local  maximum  for  the  cluster  size  of  8  observed  in 
the  computations,  whereas  in  the  experiment  it  is  ob¬ 
served  at  n  =  6.  The  reason  for  such  a  local  maximum  is 
believed  to  be  the  relation  between  the  evaporation  en¬ 
ergies  and  heat  capacities  for  the  corresponding  cluster 
sizes  (see  Table  III).  It  is  clear  that  the  values  of  these 
parameters,  used  in  the  numerical  modeling,  need  fur¬ 
ther  refinement,  and  accurate  theoretical  data  for  them 
are  indispensable  for  obtaining  better  agreement  with  the 
data.  The  agreement  between  numerical  and  experimen¬ 
tal  results  deteriorates  with  pressure,  which  is  primarily 
related  to  inaccurate  values  of  the  evaporation  energies 
and  heat  capacities  of  larger  clusters.  Note  that  no  ac¬ 
curate  quantum  mechanical  values  for  these  properties 
were  found  in  the  literature  for  clusters  larger  than  10- 
mers,  and  the  availability  of  this  information  is  critical 
for  the  present  model  (and  any  other  kinetic  model  of 
condensation) .  The  rough  agreement  of  the  data  with  an 
exponential  profile  indicates  that  the  cluster  formation  is 
mainly  based  on  monomer  addition. 


VI.  CONCLUSIONS 

A  non-equilibrium  condensation  model  applicable  for 
the  DSMC  method  is  constructed.  An  important  and 
distinguishing  aspect  of  the  present  model  is  that  all 
condensation-related  processes  are  defined  at  the  micro¬ 
scopic  level,  and  all  collision,  nucleation,  and  evapora¬ 
tion  events  depend  on  instantaneous  energies  of  colliding 
partners,  and  not  cell  temperature  or  other  macroscopic 
quantities.  The  model  was  validated  through  comparison 
with  available  theoretical  and  experimental  data  on  con¬ 
densation  rates  in  a  thermal  bath,  dimer  mole  fractions  in 
orifice  expansions,  and  cluster  size  distributions  in  nozzle 
flows.  The  model  is  based  on  a  DSMC  model  of  recombi¬ 
nation  reaction  for  the  collision  based  dimer  formation, 
and  the  RRK  model  for  the  cluster  evaporation.  Cluster 
growth  is  modeled  through  cluster-monomer  and  cluster- 
cluster  collisions.  The  energy  transfer  in  these  collisions 
is  calculated  using  the  extended  Larsen-Borgnakke  prin¬ 
ciple.  The  nucleation  rate  was  found  highly  sensitive  to 
the  values  of  heat  capacity  of  clusters. 

Two  gases  are  considered  in  this  work,  argon  and  wa¬ 
ter.  For  the  thermal  bath  relaxation,  the  present  model 
was  found  to  capture  the  dimer  equilibrium  constants 
for  water  and  argon  and  dimer  nucleation  rates  for  ar¬ 
gon  fairly  well  in  the  considered  range  of  temperatures 
from  100  K  to  350  K.  The  computed  slope  was  somewhat 
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FIG.  6.  Impact  of  the  condensation  on  the  water  translational  temperature  (K),  left,  and  axial  velocity  (m/s),  right. 
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FIG.  7.  Impact  of  the  stagnation  pressure  on  the  water  translational  temperature  (K),  left,  and  cluster  mole  fraction,  right. 


smaller  than  theoretical  for  water  clusters,  which  may  be 
attributed  to  the  accuracy  of  used  physical  and  numerical 
parameters  of  the  condensation  model. 

Comparison  of  terminal  mole  fractions  in  a  sonic  ori¬ 
fice  expansion  of  argon  with  available  experimental  data 
showed  that  the  new  model  agrees  reasonably  with  the 
measurements.  For  water,  agreement  is  reasonable  only 
for  stagnation  pressures  of  about  1.5  bar,  while  for  higher 
pressures  the  present  model  underpredicts  the  population 
of  clusters  larger  than  10-mers.  More  accurate  values  for 
the  heat  capacity  and  evaporation  energy  of  such  clusters 
are  needed  to  obtain  better  agreement  with  the  data. 
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